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Abstract 



We examine the evolution of a bistable reaction in a one-dimensional stretching flow, 



> 

in 

\o 
o 

^^ ■ as a model for chaotic advection. We derive two reduced systems of ordinary differential 

^D \ equations (ODE's) for the dynamics of the governing advection-reaction-diffusion partial 

differential equation (PDE), for pulse- like and for plateau- like solutions, based on a non- 
^ ■ perturbative approach. This reduction allows us to study the dynamics in two cases: 

first, close to a saddle-node bifurcation at which a pair of nontrivial steady states are 
born as the dimensionless reaction rate (Damkohler number) is increased, and, second, for 
large Damkohler number, far away from the bifurcation. The main aim is to investigate 
the initial- value problem and to determine when an initial condition subject to chaotic 
stirring will decay to zero and when it will give rise to a nonzero final state. Comparisons 
with full PDE simulations show that the reduced pulse model accurately predicts the 
threshold amplitude for a pulse initial condition to give rise to a nontrivial final steady 
state, and that the reduced plateau model gives an accurate picture of the dynamics of 
the system at large Damkohler number. 



PACS: 82.20.-w; 82.40.Bj 
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1 Introduction 

In natural and industrial environments active processes such as chemical or biological reactions 
are often embedded in a fluid flow. Examples include mixing of reactants within continuously 
fed or batch reactors [1,2], the development of plankton blooms and occurrence of plankton 
patchiness [3-6], increased depletion of ozone caused by chlorine filaments [7] and flame fila- 
mental structures in combustion systems [8]. The fluid flow is very often time-dependent and 
stirring. This gives rise to interesting effects, and the presence of compression and stirring in 
chaotic flows typically leads to filamentation of the chemical or biological components. This 
evidently changes the reaction dynamics. The filamentation caused by the chaotic advection 
leads to an increased surface area of the reacting components and an increase of the reaction 
output. However, if the stirring rate is too strong, the filaments become thinner and thinner 
and the reaction may stop. The corresponding flow-mediated saddle-node bifurcation has been 
reported in [9-16]. 

A natural question in such stirred reactions is: Given a chaotic flow environment and an initial 
distribution of chemical or biological reactants, is it possible to predict whether the reaction will 
take place and develop or whether it will die out? This question is of paramount importance in 
an environmental context where the reactants may be ozone in the atmosphere or pollutants 
in the ocean, or in an industrial context where one is interested in maximizing the reaction 
output. This problem is well-studied in the non-stirred case [17]. However, so far the complex 
nature of chaotic flows has prohibited an analytical treatment of this problem in the stirred 
case. 

In principle, these phenomena can be studied directly in a two- or three-dimensional reaction- 
advection-diffusion system, albeit with huge computational effort. An analytical treatment 
of the full system is prohibited by the complicated nature of the underlying equations, which 
involve multiple-scale processes. Simplified models are needed to capture essential features of 
the influence of the stirring on the reaction kinetics. In filamental or lamellar models [18], the 
two-dimensional problem of reacting tracers is replaced by a one-dimensional problem of the 

form 

d d d"^ 

—Ui - \x—Ui = Di—Ui + Ti{ui\ ki) , i = l,...,n, (1) 

for n reacting tracers {ui, i = 1, . . . ,n) with diffusion coefficients Di, reaction rates ki and stir- 
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ring rate A. Such models have been apphed by several groups, to autocatalytic, bistable and 
excitable media in several physical, chemical and biological contexts [4,5,9-14]. Phenomenolog- 
ical filamental models such as (Q) can be justified by the following consideration: The chaotic 
advection causes filaments to be stretched in one direction and compressed in another. In the 
stretched direction, the concentration is homogenized and gradients along the filaments can be 
neglected. This motivates a one-dimensional reduction for the concentration in the direction 
transverse to the filament, subject to the effect of stirring and compression. The parameter 
A can be thought of as the Lagrangian mean strain in the contracting direction, and may be 
argued to be given by the absolute value of the negative Lyapunov exponent or the (slightly 
larger) topological entropy. For a different approach to this problem see [19,20]. The validity 
of such simplified models has been numerically investigated in [21]. 

In this paper we study a bistable model with n = 1 and J-'i(ui; ki) = kiUi{ui — l){a — ui) in (jT)). 
The solutions of this model are found to be generically either of a plateau-like or a pulse-like 
character. Plateau solutions are found stably when the stirring rate is low compared to the 
reaction rate. Pulse-like solutions are found close to the saddle-node bifurcation, where the 
stirring rate is strong enough to suppress the reaction; moreover the unstable solution for weak 
stirring rates is also of a pulse-like character. The unstable solution will be of interest when 
we consider the initial-value problem and ask when an initial perturbation will develop into a 
plateau-like solution. We extend a non-perturbative approach developed in [26] to derive from 
(^ a corresponding low- dimensional system of ordinary differential equations which describes 
the time evolution of plateau-like and pulse-like solutions. We determine the equilibrium 
solutions, accurately describe the saddle-node behaviour and derive an analytical expression 
for the asymptotic width of a stationary front. By studying the phase-portrait of the reduced 
system we are able to determine the fate of various initial concentration profiles in (^; these 
predictions are then verified numerically by comparison with simulations of the full partial 
differential equation ((T)). 

In Section 121 we introduce the bistable model and discuss its solutions and their bifurcations. In 
Section 121 we describe the non-perturbative variational approach. This method will be used in 
Section m to describe pulse- like solutions and in Section El to describe plateau-like solutions. In 
SectionElwe use the reduced system derived in SectionsElandElto characterise the initial-value 
problem. The paper concludes with a discussion and an outlook for further research. 
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2 The model 

In this paper we study a bistable model 

d d d"^ 

-—u — X-—U = D-——U + Dau(M — l)(a; — n), u{x, t) ^ as a; ^ ±oo , (2) 

ot ox ox^ 

where < a < 1. The Damkohler number Da = k/\ measures the ratio of the time scales of 
fluid motion and reaction: small Damkohler numbers correspond to fast stirring/slow reaction, 
while for large Damkohler numbers the system behaves asymptotically like an unstirred system. 

The ODE corresponding to (j21) in the absence of spatial effects has two stable fixed points, 
M = and M = 1, which are separated by an unstable fixed point aX u = a. For the unstirred 
case of the PDE Q, the system is well known and well described in textbooks such as [22-25]: 
an initial perturbation which is larger than a over a finite range will spread over the whole 
domain if < a < 0.5; by contrast, if 0.5 < a < 1, an initial perturbation will decay 
to the stable state m = 0. The stirred case is much less well understood, and the initial- 
value problem has to our knowledge never been analytically treated. The stirred case was 
investigated numerically in [5,10], and semi-analytically in [15]: stationary fronts between the 
M = and u = 1 states exist as a balance between the x-dependent stirring and the diffusion- 
mediated counterpropagating fronts, for large enough values of the Damkohler number. An 
initial sufficiently large perturbation seeded at a; = spreads as a pair of fronts, driven by 
its reaction kinetics and diffusion, until the fronts reach the locations x = ±z/ where their 
velocities equal the ambient spatially dependent velocity of the chaotic stirring. 

It has been observed in [10,15] that there is a critical Damkohler number such that no stationary 
pulses exist for Da < Dae, i-e., when the time scale of the chaotic advection, Tf = 1/A, 
becomes too short in comparison with the time scale of the reaction, tr = 1/k. For large 
Damkohler numbers, an asymptotic expression for the scaling of the total concentration was 
developed in [10]. However, the techniques used in [10] cannot describe the behaviour close to 
the bifurcation point Da = Dae. Moreover, these techniques are stationary and cannot answer 
questions about the evolution of initial perturbations of reactants. 

We now describe in more detail the asymptotic form of the two nontrivial steady states for 
large Da. 
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Figure 1: The two nontrivial steady states of (J2I) for a = 0.2 and Da = 100. The 'pulse' 
solution is unstable; the 'plateau' solution is stable. 

2.1 Asymptotic steady states for large Da 



We aim to find steady states u{x) that satisfy 

xUx + Duxx + DaM(a; — u){u — 1) = 0, u{x) ^ as x ^ ±cxd, 



(3) 



when Da ^ 1. Numerical simulations show that for Da ^ 1 the stable solution is plateau-like 
and the unstable solution is bell-shaped (see Figure [T}. 

For the unstable bell-shaped solution u = V{x), we introduce 

X = Vb^x 

so that Q becomes 

DVxx + V{a-V){V-l) = 0, 

with an error of order Da^^. We may now readily determine the appropriate leading-order 
solution, which gives rise to 



V{x) = 3a{l + a + Ji(l - 2a)(2 - a) cosh 



'aDa 



x\\ + C(Da-i) 



(4) 



We note that this gives 



V{0) ~ 3a 1 + a + Ji(l - 2a)(2 - a 



-1 



(5) 
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which we shall later compare with a corresponding result obtained from a reduced test function 
ansatz using a bell-shaped test function, derived in Section HI 

In the large-Da limit, the plateau solution u = U{x) is, to leading order in Da, given by f/ = 1 
in a region —u < x < u and f/ = outside this region, where u = vDacu, for some cu = 0{1). 
The locations of the interfaces are determined by a balance between the tendency of the pulse 
to spread (with constant velocity) and the compressive effects of the imposed velocity field 
(for which the velocity is proportional to x) [10,15]; we return to this point in more detail in 
Section [5.1 .11 Around x = ±z/ there are transition regions of width (9(l/vDa). To analyse 
the region around x = i/, we introduce 



^ = vDa(a; — vDaa;). 
Then from we obtain 

D% + ujU/: + U{a -U){U -1)= 0, 
with an error of order 1/vDa. Correspondingly, 

f/(0 = i (l - tanh k/VSD] ) , (6) 



with uj = JD/2{1 — 2a). There is a similar transition region near x = —u. We shall in 
Section El use tanh profiles such as these as an ansatz for our second reduced model; since the 
ansatz captures exactly the large-Da form of the fronts, we shall find correspondingly that the 
reduced model provides an excellent approximation to the dynamics of the full PDE (P)). 

Figure Q shows the solutions of ((21) obtained by a shooting algorithm for high Damkohler 
number (Da = 100). It clearly shows that the stable solution is plateau-like and the unstable 
solution is pulse-like. Note that the stable solution for smaller Da, close to the bifurcation (not 
shown here), is also pulse-like. 



2.2 Asymptotic solutions for small Da 

We now seek an appropriate ansatz for the solution behaviour close to the saddle-node bi- 
furcation. Some insight into the form of the solution may be gained by examining the small 
Damkohler number limit of the full time-dependent problem (J2I). Beyond the saddle-node 
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Figure 2: The two nontrivial steady states of Q for a = 0.2 and Da = 10; tlie upper and lower 
solutions are, respectively, stable and unstable. Shown in each case with crosses are fits to a 
Gaussian. 

bifurcation point, i.e., for Da < Dae, no steady-state solutions exist, and any initial condition 
decays to zero. In fact, if we expand the solution as a series in the Damkohler number according 
to u{x,t) = Da,UQ{x,t) + (9(Da^), the equation is satisfied at leading order by 



Uo{x,t) = /o(t)exp(-(wo(t)x)^ 



(7) 



where wq -^ 1/vZD and /o — > as t ^ cxd; the latter limit reflects the nonexistence of 
nontrivial steady states for Da < Da^. Although Dae = C^(IO), which is clearly not small, it 
seems that the leading-order term ((7j) provides an accurate representation of the solution to 
(121) close to the saddle-node point Da = Dae - see Figure |21 So in the first instance we shall 
approximate both stable and unstable solutions close to the saddle-node point by a Gaussian; 
we shall later consider, in somewhat less detail, other possible profiles of 'bell-shaped' form. 

Using a non-perturbative method developed in [26] we now derive a set of ordinary differential 
equations which describe the dynamics of these plateau- and pulse-like solutions. In the next 
section we briefly explain this variational approach. 



3 Nonperturbative, variational method 



A method was developed in [26] to study critical wave propagation of single pulses and pulse 
trains in excitable media in one and two dimensions. It was based on the observation that close 
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to the bifurcation point the pulse shape is approximately a bell-shaped function. Numerical 
simulations and the asymptotic analysis in Section 12.11 and 12.21 show that this is also the 
case for the bistable model (J2I) close to the saddle-node bifurcation at which the stable and 
unstable solutions are born, and on the unstable branch at large Damkohler numbers. A test- 
function approximation that optimises the two free parameters of a bell-shaped function, i.e., 
its amplitude and its width, allows us to find the actual bifurcation point, Dae, and determine 
the pulse shape for close-to-critical pulses at Damkohler numbers near Dae. We note that the 
framework of asymptotic techniques, such as inner and outer expansions where the solution is 
separated into a steep narrow front and a flat plateau, are bound to fail close to the bifurcation 
point, because the pulse is clearly bell-shaped, and such a separation is no longer possible. We 
shall make explicit use of the shape of the pulse close to the critical point and parameterise 
the pulse appropriately. 

We choose u{x, t) of the general form 

u{x) = fit)U{T]) with 7] = w{t)x , (8) 

where U{ri) is a symmetric, bell-shaped function of unit width and height, and f{t) is the 
amplitude of the pulse. Close to the saddle-node bifurcation and for the unstable solution at 
large Damkohler numbers, the solution is indeed a bell-shaped pulse, which we approximate 
by a Gaussian (see Figure |21). However, our results do not depend on the specific choice of the 
test function, and the numerical values differ only marginally when sech-functions are used 
instead, for instance (see Figure Ej). We restrict the solutions to a subspace of a bell-shaped 
function /[/(r/), which is parameterised by the amphtude, fit), and the inverse pulse width, 
w{t). The evolution of these parameters is then determined by minimizing the error made by 
the restriction to the subspace defined by (jS)); this is achieved by projecting onto the tangent 
space of the restricted subspace, which is spanned by du/df = U and du/dw = frjUrj/w. We 
set to zero the integral of the product of Q with the basis functions of the tangent space (over 
the entire ?7-domain). This leads to ordinary differential equations for f{t) and w{t)] it also 
yields an approximation to the critical Damkohler number Dae. 

One may in fact choose any plausible test function and optimize the restriction to the particular 
space of solutions by varying the parameters. This variational idea is not restricted to pulse- 
like solutions. Far away from the saddle-node bifurcation, at large Damkohler numbers, where 
the stable asymptotic solution comprises a pair of fronts (see Section EH} we may use as a test 
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function a combination of tanh-functions 

u{x, t) = |/(i) [tanh(r7 + a{t)) — ta.nh{ri — a{t))] with r] = w{t)x, (9) 

where a{t) = w{t)v{t) and the fronts are located around x = ±z/. In a manner analogous 
to that described above for a pulse ansatz, the variational approach allows us in this case to 
determine the time evolution of the solution amplitude f{t), the inverse interface width w{t) 
and the front locations ±z^(t) by projecting Q onto the tangent space of the restricted solution 
subspace, which is now spanned by du/df, du/dw and dujdv. 



4 Stable and unstable pulse solutions 

Near the turning point at Da = Dac(a), both nontrivial steady-state solutions take the form 
of bell-shaped pulses near x = 0. We note that in the literature a •pulse commonly refers 
to a stable homoclinic solution; however, here we use the term merely to mean a bell-shaped 
function, regardless of whether it is stable or unstable. This terminology allows us to distinguish 
between the two test function ansatze (jH)) and Q. 

To determine a reduced model for such pulse solutions, we therefore introduce the ansatz (jH} 
and to determine the evolution of /(t) and wii) we enforce the vanishing of the two inner 
products 

((■Uj — XM^ — DWa-a; — Da «(« — «)(«— 1))-Uj) = 

(^{ut — xu^ — Duxx — ^^u{a — u){u—V))u^ = 0, 



where 






For a Gaussian pulse, with U{ri) = e '' , we note the requisite integrals 

<-)-©"■ <-"^>-(i)"' <^^-'^)^i©"- 

Thus after some algebra we find that f{t) and w{t) evolve according to 

/ = /|-Daa-2DM;^ + |Da(l + a)yi/-|Day^/2| (10) 

jl-2Dw;^ + iDa(l + a)y|/-|Day|/2|. (11) 



w = w 
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In order to examine the dynamics of this system in the phase plane it is useful for brevity of 
notation to write it in the form 

/ = f {-1^1- l^2W^ + fisf - l^if} (12) 

w = w{X^-\2W^ + \sf-\,f}, (13) 

and note that the /i, and Aj are all positive; also fj,2 = A2 (this equality does not necessarily 
hold for other profiles U{ri)). 

We now examine the steady states {w,f) = {w, f) of (fT2|) . |T3|) and their stability. Note 
that the steady-state problem has been considered previously [15], although only the solutions 
identified below as 'P4' were considered there. Although not all steady states of ()12j). ()13|) 
correspond to distinct physical solutions, they help us to map out the structure of the phase 
plane, and hence determine the dynamics of (fT2|l . (fT3|) . 

Pi: (17,7) = (0,0) 

This steady state exists for all parameter values, and a linearisation about Pi yields the eigen- 
values Ai and — /ii. Therefore Pi is a saddle point, unstable to perturbations in w and stable 
to perturbations in /. 

P,: (^,7) = ((Ai/A2)i/^0) 

This steady state exists for all parameter values. A linearisation about this point yields the 
eigenvalues — /xi — ^12^1/^2 < and — 2Ai < 0. Thus P2 is a stable node. 

P3: 117 = 0,77^0 



Here 



/^4/ -/U3/ + AH = 0. (14) 



This quadratic has two real, positive solutions 

J ^ /i3 - (/ij - 4/i4/il)^/^ 
2//4 
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for all parameter values; indeed /_ and /_,. are independent of Da. Note that, since P3 corre- 
sponds to infinitely wide pulses (W = 0), the roots of the quadratic would be a and 1 if there 
were no approximation due to the Gaussian pulse ansatz. A linearisation about either steady 
state yields the eigenvalues 

ei = (/is - 2/14/)/, 62 = (Ai + As/ - A47 ). 

Now /is— 2/i4/_|_ = =F(/i|— 4/i4/ii)^/^, so the eigenvalue ci indicates stability for /+ and instability 
for /_ (ei corresponds to perturbations in /, with w = 0). The eigenvalue 62 corresponds to 
disturbances in w and satisfies 

/i4e2 = (A3/i4 - A4/i3)/ + (A4/ii + XijJ'a); 

it seems to indicate instability for /+ and either stability or instability for /_, depending on 
the parameter values. 

P4: w^O,J^O 

Here / satisfies 

(/i4-A4)7'-(/xs-As)7 + /xi + Ai = 0. (15) 

We note the coefficients (cf. [15]) 

3Da ^ 5Da(l + a) ^ , , 

fi4- M = —y=, /is - As = —y= , /ii + Ai = 1 + Daa. (16) 

The quadratic (|T3|) has two real, positive roots provided 



a <am = ^ 0.47448, where q - 



2q • ' 81v^' 

and 



Da > DaP'^^'^" = \q{l + af - a 



-1 



Note that am is slightly smaller than 0.5, the threshold for existence of non-zero solutions in 
the unstirred case. 
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In this case, we denote the two solutions of (fT3j) by /+>/*. Correspondingly W = w^, where 

/i2W±^ = -/i4/±^ + /is/i - /ii- (17) 

Since the right-hand side of ()17|) must be non-negative, we deduce, by comparison with (fT^ . 
that each of /^ and /^ lies between /_ and /+. At Da = Da^?"^*^, the two solutions collide in a 
saddle-node bifurcation; no such solutions exist for Da < Da^*^^*^*^. 

A linearisation about either of the P4 steady states yields, for / = f + Sf(t) and w = w + 6w{t), 
that 

Sf = (a^s/ - 2/14/ ) 6f - 2/i2w7 6w 

Hence the growth rate a satisfies 

a^ + ba + c = 0, (18) 

where 

b = 2A2U7^ - fij + 2fxJ = -2/ii + /ig/ 

and 

C = 2/i2W^7(2(/i4 - A4)/ - (/is - A3)). 

Note that for /^, c > 0; for /^, c < 0, indicating instability. The bifurcation point c = 
corresponds to (the limits uJ = or / = or) Da = DaJ?" '^'^. Note that a = (9 (Da) as Da -^ 00. 

Two phase planes for (J12p . (jl3|) illustrating the saddle-node bifurcation at Dae are shown 
in Figure El If Da < DaJ^'^^^, the fixed points P4 do not exist, and all initial conditions off 
the /-axis asymptote to the solution {w,0) at large time (physically the trivial state u = 0). 
If instead Da > DaJ?"^*^, the P4 solutions exist, and there is a separatrix that divides initial 
conditions leading to the physically trivial state {w, 0) from those leading to the physically 
nontrivial state {w*^_, f^). 

Figure ID shows a comparison between results from the Gaussian pulse model and steady-state 
solutions of (J2) computed using shooting. The variational result (fTS|) accurately reproduces 
the bifurcation behaviour, both in terms of the critical Damkohler number and the generic 
quadratic behaviour close to the saddle-node bifurcation. A motivation for our Gaussian 
ansatz was given by the small-Da asymptotic argument in Section 12.21 However, in Figure El 
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separatrix 




Figure 3: Two of the possible phase planes for (fT^ . (fT!^ illustrating the saddle-node bifurca- 
tion. Left: Da < Da^^''^ and the fixed points P4, do not exist. Right: Da > DaP"^'^ 



X! 




Figure 4: Comparison between the steady states of Q computed by shooting (solid line), the 
pulse model (fT3|) (dashed line) and the plateau model (IHU)) (dotted line). Here a = 0.2. 

we show that good agreement is also given for different ansatze such as U{t]) = sech"?;. We 
thus conclude that the actual form of the pulse shape is not particularly important (and hence 
we cannot expect to deduce this pulse shape from any asymptotic theory). 



By contrast, the asymptotic approaches which are usually employed, and which assume plateau- 
like solutions with well-separated constant plateaus, fail to describe the bifurcation behaviour. 
We note that the agreement between the pulse model and the full problem (J2I) is excellent on 
the lower branch. Figure IHl shows the excellent agreement between m(0) = raaxu{x) computed 
according to the pulse result (fT3|) in the limit Da -^ 00 and a direct large-Da asymptotic 
analysis of (0), given by (0)- 
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Figure 5: Comparison between different shaped pulses in the model p5|l . The solid curve 
represents the steady states of © computed by shooting; the dashed lines represent solutions 
of the pulse model ()15|) for test functions U{ri) = e~^ , sech^?7 and sedx^r] (respectively from 
right to left at the saddle-node bifurcation point). Parameters as in Figure 01 




Figure 6: Asymptotic value of maxM(x) for the unstable steady state at large Da. Solid line: 
result from (jlSp . Dashed line: result from (jSJ. 
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5 Plateau solutions 

Away from the saddle-node bifurcation, the stable solution takes the form of a nearly uniform 
region around x = 0, surrounded by a pair of fronts. This solution may be captured by writing 
u{x,t) in the form 

u{x,t) = \f{t)(p{r]) with r] = w{t)x, (19) 

where 

0(77) = tanh(?7 + a{t)) — tanh(?7 — a{t)). 

Then the evolution of /(t), w{t) and z/(t) = a/w (or, equivalently, /, w and a) may be 
determined by forcing the inner products of (j2)) with each of uj, u^ and u^ to be zero, where 

du 



«/ 



u„ 



M,. 



df 
du 



Uiv) 



u 
1 



■TT- = \f (2; + i^) sech^(?7 + a) — ix — u) sech^(r7 — a) 
ow ^ L 



du 
du 



with 



kfw-^iv), 



ip{ri) = sech (77 + a) + sech (r] — a). 



The corresponding steady-state calculation has been considered in [15]; here we extend their 
analysis to the time-dependent problem. This permits us to determine which initial conditions 
lead to excitation and which to extinction. Details of our derivation, including analytical 
expressions for all the requisite inner products, are given in the Appendix: we find that the 
equations for the time evolution of /, w and u take the form 



Xiwf + Aa/w + Xsw'^fu 
(3iwf + (32fw + (3^w^fu 



wf [A4 + X^Dw' + Da(A6 + Xjf + Xsf) 
wf /i4 + /is-Dw^ + Da(/i6 + firf + fisf] 
wf [/?4 + f3,Dw^ + Da(/56 + /?7/ + (38 f) 



(20) 
(21) 
(22) 



The various coefficients Aj, /Xj and /3j are complicated functions of a and a, and are listed in 
the Appendix. 
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5.1 Steady states 

We shall not attempt here a full analysis of the three-dimensional phase space for /, w and v\ 
instead we consider only nontrivial steady-state solutions for which /, w and v are all nonzero. 
We expect two such solutions when Da exceeds some threshold, one stable and one unstable. 
We shall consider the dynamics of ()2()j) - (P^ later. 

The equations for the (nontrivial) steady states may be written in the form 

A4 + As/^w^ + Da(A6 + A7/ + As/') = (23) 

/i4 + Ais^w' + Da(/i6 + /^t/ + /is/') = (24) 

(3^ + P^Dw^ + T)^{(3<, + (3jf + (3^P) = 0. (25) 

These equations are readily solved without recourse to nonlinear root finding, for example in 
the following way. First we fix values for a and a. Then we use fj23|) to write 

2 -A4 - Da(A6 + A7/ + As/') .^„. 

w = ^ . (26) 

Then, upon substitution of this expression for w"^ in ^^, we have 

/i; + Da(/i; + /i'7/ + /i^/') = 0, where fi[ = ^,-^K (27) 

Likewise 

/?; + Da(/3; + /?;/ + /?;/') = 0, where p[ = p,-(^X,. (28) 

Now by eliminating Da between (J2ZI) and (|2H1), we find the following equation for /: 

(/il/5^ - ^',P',)f + {ix',P', - ^',P',)f + fi',P', - K/3; = 0, (29) 

which is of the form 

A{a)f + {l + a)B{a)f + aC{a) = 0. (30) 

The quadratic equation ()30|) needs careful interpretation. First we note that the coefficients 
A, B and C are functions of a alone, and hence for given values of a and a there are in general 
two solutions of (|29|) or none. When there are two solutions, they correspond, through (|27p or 
()28|). to distinct values of Da; numerically, we observe that one value of Da is positive and the 
other negative (and hence physically unreasonable). Furthermore, when we reconstruct from 
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Figure 7: Steady-state solutions of the plateau model (j^ - lj^ (solid line) for a = 0.2. The 
dashed lines give large-Da asymptotics (|H^. (jHHj) of the plateau model for the upper branch. 
The asymptotic results are indistinguishable from the numerical results for / and a. 

(j26|) the corresponding values of w, we find numerically that the root of (j29j) corresponding to 
positive Da seems always to give rise to w'^ > 0, and so is indeed physically reasonable. Thus 
we emphasise that both branches of solutions evident in Figures [7| and |H1 arise from one of the 
roots of (J30|) for fixed a, as a is increased monotonically. 

The steady states for a = 0.2 (see [15]) are shown in Figure [7| Note that solutions exist only 
for Da > Da^'^*°^"(a). Figure |H1 shows steady-state solution branches for a range of values of 
a. The excellent agreement between solutions of the plateau model and steady states of (j21) on 
the upper branch is illustrated in Figures E] and IHl Of course, we expect a tanh-test function 
ansatz to fail close to the bifurcation point Dae where the solution rather is bell-shaped and is 
accurately described by a Gaussian ansatz as in Section EJ 

Of particular interest is the limit of large Damkohler number, Da ^ 1. The numerics suggest 
the following scalings. On the upper branch of solutions 

a = aDa, w = wDa^'^, / = /, (31) 

where a,w, f = 0{1); on the lower branch 

a = a, w = wDa^'"^, f = f, 

where d,w, f = 0{1). 
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Figure 8: Steady-state solutions of the plateau model (pHjl - fOK|l for a = 0.1,0.2,0.3,0.4. (The 
corresponding value of 10a is indicated next to each curve.) 

5.1.1 Large Damkohler number, Da ^ 1, upper branch 



We now construct the plateau solutions on the upper branch, in the limit Da ^1. In that 
limit, we find 



A4 ~ —a 




/^4 -Ts^^-l 


A5~-| 




/i5~-| 


Ag ~ —2aa 




/i6 ~ -|a 


A7~2(l + 


a)a 


/i7~ i(l + a) 


Ag ~ —2a 




/^8~-i 



(34 


r\j 


-¥ 




/?5 


^ 







/36 


r\j 


—a 




/?7 


r^ 


1(1 + 


a 


Ps 


rvj 


1 

2' 





Recall that a = (9(Da). The result for /?5 indicates that this quantity is exponentially small in 
Da. Thus at leading order in Da we have, from ()23 p - (P3j) . 



-a+{l + a)f-p 



11 -f2\ 



-iDu;' + i(-a + (l + a)/-ii/ 



-fa + (-a + |(l + a)/-i/2) 





0. 



The first of these equations gives / = 1 or / = a; the latter solution must be rejected since the 
last equation would then give a = —{2 — a)a/4 < 0. The choice / = 1 is consistent with our 
numerical calculations, which indicate that this is the solution relevant to the upper branch. 
With / = 1, the remaining equations give Dw"^ = | and a = 1(1 — 2a). Thus in the large-Da 
limit. 



/ 



w 



Da\i/2 
8D. 



a~i(l 



2a)Da. 



(32) 
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Figure 9: Comparison between stable steady states of the plateau model (crosses) and of 
the advection-reaction-diffusion equation (j2)) (solid line). Here a = 0.4 and Da = 100.0359 
(implied by our choice of a); for the plateau-model solution, a = 5.0, / = 1.000179 and 
w = 3.62460. 



Note, as a consequence, that for large Da we have 

V2DDa(i-a) 



u ~ 



(33) 



in agreement with the phenomenological argument given in [10, 15]. This formula can be 
understood if we note that the location of a stationary front is given approximately through a 



balance between the front velocity vq = \^2DDa (| — a) of the unstirred problem [10] and the 
velocity of the chaotic stirring (i.e., x in the scaling of Q). Thus in the stirred problem, the 
front has zero velocity when vq = x, which implies Vq = u, and (j33|) follows. 



Finally, the solution takes the form 



M ~ - < tanh 



'Da (l-2a)Da 

Sd""^ 4 



— tanh 



'Da (l-2a)Da 
Sd"" 4 



(34) 



A comparison between this expression and a numerical calculation of the stable steady state 
of (0) is shown in Figure 01 to graphical resolution, the two graphs are indistinguishable. 

Figure UHl shows how well the ansatz (jTHI) captures the width of the steady-state solution. 
Plotted is 2(t'approx — t'exact), whcrc z/approx and Inexact a^G, respectively, the half-widths of the 
solutions to the plateau model and the full problem (0). In each case only the stable solution 
on the upper branch is considered, and the 'half-width' is defined as z/approx, t'exact = a;+ > 0, 
where u{x-^) = 0.5. The difference 2 ( z/approx— z^cxact) clearly decreases with increasing Damkohler 
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Figure 10: Comparison between the width of the steady-state solution as obtained from the 
plateau model and from the advection-reaction-diffusion equation (j21). Here a = 0.4. Shown 



is 2(iy. 



approx 



t'cxact) against the Damkohler number Da. 



number, which is expected, because the tanh-front ansatz gets closer to the actual solution, 
which asymptotically is a pair of tanh- fronts (see Section l?!T|) . 

Solutions of the plateau model on the lower branch are of less interest because they do not 
correspond closely to steady-state solutions of the full PDE, where the unstable solution is 
pulse- like. 



5.2 Stability 



Suppose that {f,w,v) is a nontrivial steady state of (j2ni)^(E2I) (i-e., f.w.i' ^ 0). Since the 
various coefficients in (j^Ujl - lj^^ are functions of a, it is useful to consider the linearised evolution 
of 5f, 5w and 5a, where in the linearisation 5a = w5v + v5w. Thus we have 

\iw5f + \2f5w + \zfw{5a - v5w) = Da fw{X7 + 2Xsf)5f + 2Dfw'^\^5w 

+ fw [OXa + Dw'^d\ + Da((9A6 + fdXj + fdX^)]5a 



where 



dX,= 



da 
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and where all nonperturbation quantities are evaluated at the fixed point. There are two 
further equations, with Aj t— > fii and Aj h-> /3j. The three may be written together as 



A 



6w 
\ 6d / 



B 



6w 

V ^« / 



where the constant-coefficient matrices A and B are readily determined from the three pertur- 
bation evolution equations. 



5.2.1 Large Damkohler number, Da ^ 1, upper branch 



On the upper solution branch, in the limit Da ^ 1, we have 



/2 



A 



aw 



—a w 



\ 



\ w 



\w W'-Q) 
-fa Iw ) 



/-2Dawa(l-a) -\Dw'^ -w\ 



B 



^Daw(6Q;-5) -\Dw^ 



\ -|Daw(l-2a) 







\w j 



where the zero entries indicate terms that are exponentially small in Da. Thus, in view of the 

scalings (jSH), 

/ 2a«;Da^/2 -aDa wV)^l^ \ 



A 



i«;Dai/2 ^(7r2-6) 







V wDa^/2 -faDa \wV)^l'^ j 



and 



/ -2wa(l - a)Da^/2 -\Dw^X)& -wDa^/^ \ 



B 



^«;(6« - 5)Da^/2 -Ww'^V) 



\ -|w;(l - 2a)Da^ 



-|wDa^/2 ) 



It then follows (after some algebra) that 

/ 



A-^B 



-(l-a)Da 
Da^/2 



2(7r2-6) 



3a 



V 2(7r2-6) 



Da^ 



XIDwa 0/9 
-Da^/^ 



TT^ — 6 



/ 



We require the dominant eigenmodes of A ^B. There are three eigenvalues: 

Ai 1 and A2,A3 = C(Da). 



(35) 
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The eigenvector corresponding to Ai has Sf = Sw = and Sa = 1, say; thus disturbances in a 
decay hke e~*. We find A2 ~ —(1 — a;)Da; recaUing the scahngs / = 0{1), w = (9(A/Da) and 
u = 0(-\/Da), we find that the associated eigenvector has {Sf,Sw,6a) ~ {F,WDa}''^, ADa), 
with 

The final eigenvalue is 

3 

A3 ~ -— rDa, 

^ 2(7r2 - 6) ' 

with eigenvector {6f,6w,6a) ~ {F,WDa!^''^, ADa^), where 

^ 1 - a I F = -'^^W and A = SDwaW. 



\ 2(7r2-6) y 3a 

Note that this calculation confirms the stability of the upper branch solution, since Ai^2,3 < 0. 
Furthermore, the scalings (jH^j) confirm the observation from numerical simulations of the PDE 
(12)) that the relaxation of the front separation to its equilibrium value is slower than the 
relaxation of either the amplitude of the solution or the width of each front. (Provided the 
fronts are reasonably well separated initially, numerical integrations of either the PDE or the 
system ((201)^(1221) indicate that first / adjusts to approximately its long-time value; then w 
adjusts to give each front the appropriate width, and finally a adjusts so that the fronts are 
separated by the correct distance.) 

Since < a < | and since |(7r^ — 6)~^ ~ 0.3876, we have 

IA2I > IA3I > |Ai|. 



6 Comparison between dynamics of the reduced models 
and the full PDE 

Our models extend the steady-state analysis of [15] by providing approximations to the dy- 
namics of the PDE (j21). This allows us to examine the evolution of initial conditions which 
are either plateau-like or pulse-like. For instance, using the Gaussian pulse ansatz, the dy- 
namics is approximately reduced to that of the two-dimensional system (fT^ - (fT!?|) . Given an 
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/ 




Figure 11: Comparison between the separatrix / = ff^^^'^{w) determined from the Gaussian 
pulse ansatz (continuous hne) and the threshold / = ff^^{w) obtained from simulations of 
the full PDE Q (crosses) with initial profile u{x, 0) = fe~^^^^ as a function of w; parameter 
values are a = 0.2 and Da = 16. 

initial pulse-like solution we may now solve the initial-value problem at Damkohler numbers 
close to the critical Damkohler number Dae, where both stable and unstable solutions are well 
approximated by a bell-shaped function. We can readily compute the corresponding separa- 
trix in {w, /)-space between those solutions that asymptote at large times to the uniform zero 
solution, and those that asymptote to a nontrivial steady state. We denote the separatrix by 
/ = /f"'^°(w^)- For initial conditions (w(0),/(0)) below the separatrix the solution will die 
out, whereas for initial conditions above the separatrix the reaction will continue and develop 
into a stable stationary pulse. It is of particular interest to determine the extent to which 
this result provides useful corresponding information for the full PDE (0). We have therefore 
carried out a large number of simulations of the PDE, each from a Gaussian initial condition 
u{x,0) = fe~^^^^ , in an effort to compute a corresponding separatrix: / = ff^^{w). For a 
given value of w, we use bisection to determine approximately the corresponding value of / that 
delineates the two large-time behaviours. The results, and a comparison with the separatrix 
from the Gaussian pulse ansatz, are given in Figure ITTl The agreement is excellent. There is 
similar agreement at higher Da; for example, at Da = 100, and w{0) = 1, the Gaussian ansatz 
gives a threshold /(O) = 0.2537, while from the PDE we find 0.2525 ± 0.0005. 



Assessing the usefulness of the plateau model is rather harder. This is because if we start with 
an initial profile of the form (fTIUl . with well separated fronts, and u{x, 0) approximately uniform 
between those fronts, then in the uniform region the initial evolution of u is well approximated 
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Figure 12: Comparison between numerical simulations of the PDE Q and the plateau model 
(|2(jp - (j22j) . Top: initial and final states, with arrow indicating direction of temporal evolution. 
Bottom left: u{0,t) against time. Bottom right: 'width' of the pulse, 2z/. In each of the lower 
plots, PDE results are given by a solid line and plateau-model results by crosses. 



by the ODE 



du 

dt 



Da.u{a — u){u — 1) 



so that / is initially decreasing if /(O) < a and increasing if /(O) > a. Thus the threshold 
value of /(O) leading to zero or nonzero large-time states is rather trivially /(O) ~ a. 

The agreement between numerical integrations of the PDE (0) and the ODE system for a 
plateau-like solution ()20 |) -()22 |) is illustrated in Figure ^1 which shows simulations from an 
initial plateau-like state of the form (jl9|) . for a = 0.2 and Da = 100. The initial condition 
is / = 0.22, w = 1 and a = 2.5. Two quantities are plotted: the maximum u{0,t) and the 
'width' of the solution, 2i/, defined again as 2x+, where u{x+,t) = 0.5 and x+ > 0. (Note that 
initially the entire solution lies below 0.5, so the 'width' is initially undefined.) The solid lines 
represent numerical simulations of the PDE (j2)); the crosses represent simulations of the ODE 
system ()20|) - ()22|) . The agreement is excellent, indicating that at large Da the plateau model 
accurately describes not only the stable steady state, but also the dynamics of solutions to the 
PDE. 
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The agreement between the plateau model and the full PDE is remarkable because, as shown 
in Section im the unstable solution is pulse-like and not plateau-like. Furthermore, as shown 
in Figure m the unstable solution predicted by the plateau model is quantitatively incorrect. 
However, some insight into the agreement may be obtained by noting that the unstable pulse- 
like steady-state has eigenvalues a = (9(Da), which are large in the limit Da —>■ oo, and which 
indicate that the associated instability is reaction-dominated. Our simulations of the PDE Q 
confirm that the initial condition illustrated in Figure IT^ first rapidly undergoes an adjustment 
of the central concentration to m ~ 1, consistent with a time scale (9(1/Da), then a slower 
relaxation of the front width and location. 



7 Conclusions 

We have developed reduced systems of ordinary differential equations describing the time 
evolution of both plateau-like and pulse-like solutions of a one-dimensional bistable reaction- 
diffusion partial differential equation for a stirred flow environment. The PDE possesses a 
trivial steady state in which all the reactant is used up, and, when the dimensionless reaction 
rate Da is great enough, a pair of nontrivial steady states, one stable and one unstable, born 
in a saddle-node bifurcation at Da = Da^. The pulse model was shown to provide a good 
approximation close to the bifurcation point and for the unstable steady solution at large Da. 
The plateau model applies for large Da, and captures well the form of the stable steady state. 
By studying equilibrium pulse solutions, we were able to accurately describe the bifurcation 
behaviour of the system near Dae. For large Damkohler numbers our reduced plateau model 
allowed us to describe the stationary fronts, and derive a commonly used phenomenological 
formula for the width of a stationary front. 

The time-dependent ODE models allowed us to study the initial- value problem for Q. We 
found that although the unstable solution is actually of a pulse-like shape, the dynamics of 
(j2)) is well captured by the plateau model for high Damkohler numbers. This remarkable and 
nonintuitive result is due to the fact that the growth rate of the instability of the unstable 
solution is (9 (Da), whereas the fastest growth rate for the shape-change of a front is 0{1). Our 
numerical solutions of Q seem rapidly to approach the plateau form, and thus the dynamics 
of the PDE (j2)) is described accurately by the plateau model. 
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Besides being instructive in themselves, one-dimensional models such as (0) aim to provide 
information about the evolution of chemical reactions in a two- or three-dimensional chaotic 
flow. It is now known [21] that simple variants of reduced models perform poorly in predicting 
the yield of multi-stage reactions, and so a more sophisticated analysis, as presented here for 
the bistable case, is therefore warranted. One particularly promising extension [27] of one- 
dimensional models such as (P) involves two stages. In the first stage, a point x in the two- 
or three-dimensional problem and a time t are selected, at which the chemical composition 
is required. This point and its associated stable manifold are then advected backwards in 
time, to t = 0, and the initial condition for the system is imagined to be sampled by the 
stable manifold (which is now exponentially stretched and highly contorted). The second stage 
involves carrying out a one-dimensional simulation of the advection, reaction and diffusion 
along the evolving stable manifold (now forwards in time), using a model such as ([Q), but with 
a nonconstant compression rate 
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Appendix 



Here we record some details of the derivation of the coefficients in the plateau model ()2Up^ 
). We begin by noting that wuw = jfri(j)'{ri) + ucp^ = ^/{rjcj)' + wuip), xUx = ^fv^P'iv) ^^^ 



u. 



1 fn„'2All 



fw^"{r^). 



Now we list the various inner products that are needed for the calculation: 



{UuUf) 

{<f>"uf) 
{u{a — u){u — l)uf) 

{<) 
{u{a — u){u — 1)mui) 
{u(a — u){u — V)Ui,) 



1/^2 
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4w 

lfw{ip(f)) 



8m; * 



^') + ^(^0) 
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-f 
f 



{rfct)') + 2wz/(r/0V) + w^v^{'il)^) 



i/'(#V) + \fwv{i,') 



.2 ,/2\ 



/^, 






d 



Aw \ da 



(,'2\ 



1 



a(02) + i(l + «)/(0')-T^/'(0') 



yv{^') 



\fw{ri(t)'i)) 



Ifw^ [-|«(0') + 1(1 + «)/(0') - Tefi^') 



All of the requisite inner products can be found exactly: 

{(f)'^) = 4(2acoth2a-l) 

1 d 

(^0) = _ (02^ ^ ^ ^^^Yi_ 2a - 8a cosech^2a 



2 da 



j,'2 



'1 + 3 cosech 2a — 6a cosech 2a coth 2a) 



') = 4(4a coth^ 2a + 2a cosech^2a - 3 coth 2a) 
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') = f (-11 + 12acoth^2a - 15cosech^2a + 18acoth2acosech^2a) 






^1 — 3 cosech^2a + 6a cosech^2a coth 2a) 



1^2 2 I 4^2\ 4:^/2 



2(i7r^- f + |a j - |a(7r^ + 4a^)coth2acosech^2a + |(7r^ + 12a^)cosech^2a. 



The various coefficients in (J2(J|) - (J22|) are then given by 



Ai 



1/^2 



) , yUi = As 



1 /A'i 



) + \a{^|J<p), (3, = X, = l{^p<p) 



l/J,/2\ 



A4 = -|(</)^), A5 = -|(r), A6 = -ia(02), A7 = |(l + a)(0=^), Ag 



(r/V') + 2a(#V) + a2(V'') , /^a = /is = ^ (#V) + a(V'') 



/^2 = 4 _ 

/^7 = -^(1 + a) U - a^ j (0') , /is = ^ f 1 - a^ j (0^) : 



rO; 1 — a- 



9a 



/?6 



l«l:('^')' /?v = ^(i + «)|-(0^), /?, = -^|-(0^). 
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da 



da 
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